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ABSTRACT 

A self-consistent formalism to jointly study cosmic reionization and thermal history of the in- 
tergalactic medium (IGM) in a ACDM cosmology is presented. The model implements most 
of the relevant physics governing these processes, such as the inhomogeneous IGM density 
distribution, three different classes of ionizing photon sources (massive PopIII stars, PopII 
stars and QSOs), and radiative feedback inhibiting star formation in low-mass galaxies. By 
constraining the model free parameters with available data on redshift evolution of Lyman- 
limit absorption systems, Gunn-Peterson and electron scattering optical depths. Near InfraRed 
Background (NIRB), and cosmic star formation history, we select a fiducial model, whose 
main predictions are: (i) Hydrogen was completely reionized at z « 15, while Hell must have 
been reionized by z « 12, allowing for the uncertainties in the ionizing photon efficiencies of 
stars. At z « 7, Helll suffered an almost complete recombination as a result of the extinction 
of PopIII stars, as required by the interpretation of the NIRB. (ii) A QSO-induced complete 
Hell reionization occurs at z = 3.5; a similar double H reionization does not take place due to 
the large number of photons with energies > 1 Ryd from PopII stars and QSOs, even after all 
PopIII stars have disappeared. ( Hi) Following reionization, the temperature of the IGM corre- 
sponding to the mean gas density. To, is boosted to 1.5 x 10'* K; following that it decreases 
with a relatively flat trend. Observations of Tq are consistent with the fact that He is singly 
ionized at z > 3.5, while they are consistent with He being doubly ionized at z < 3.5. This 
might be interpreted as a signature of (second) Hell reionization. (iv) Only 0.3% of the stars 
produced by z = 2 need to be PopIII stars in order to achieve the first hydrogen reionization. 
In addition, we get useful constraints on the ionizing photon efficiencies (which are combina- 
tion of the star-forming efficiency and the escape fraction of ionizing photons from collapsed 
haloes) of PopII and PopIII stars, namely, epopii < 0.01, 0.002 < epopiii < 0.03. Varying the 
efficiencies in these two ranges does not affect the scenario described above. Such model not 
only relieves the tension between the Gunn-Peterson optical depth and WMAP observations, 
but also accounts self-consistently for all known observational constraints. We discuss how the 
results compare with recent numerical reionization studies and other theoretical arguments. 



1 INTRODUCTION 



The study of cosmic reionization has witnessed a very strong ad- 
vance in the last few years (for a recent review see Ciardi & Ferrara 
2004), thanks both to the availability of QSOs located at redshifts 
> 6 which allow to probe the physical state of the neutral part of the 
intergalactic medium (IGM) via absorption line experiments (Fan 
et al. 2000; Fan et al. 2001; Fan et al. 2002; Fan et al. 2003), and to 
the determination of the Thomson electron scattering optical depth 
Tel from the first year WMAP data (Spergel et al. 2003). A consid- 
erable tension from these two data sets has arisen concerning the 
epoch of hydrogen reionization. In brief, while the rapid increase 
in the (Gunn-Peterson) Lya opacity at z > 5 could be interpreted, 
to a first analysis, as an indication of a reionization occurring at 
z ~ 6, the large inferred values of Td would imply a reionization 
epoch at z > 10. In this regard one has to recall that (i) the first 
year WMAP result has large error-bars; (ii) the limits on t^i de- 
pend on the exact analysis technique employed; (iii) the robustness 
of the result needs to be confirmed using more data, while on the 
other hand (i) the Gunn-Peterson test is very sensitive to even tiny 
amounts of neutral hydrogen, resulting only in a lower limit to the 



neutral hydrogen fraction, xm ^ 0.1% (White et al. 2003); (ii) the 
experiment has been carried out on a handful of QSOs above z = 5 
and therefore its statistical significance has still to be evaluated; (iii) 
the abrupt increase might be in part due an incorrect conversion of 
Ly/3 to Lya optical depths (Songaila 2004). Nevertheless, reasons 
for the apparent discrepancy must be considered seriously. 

A crucial issue about reionization is that this process is tightly 
coupled to the properties and evolution of star-forming galaxies and 
QSOs. Hence, the first requirement that any reionization model 
should fulfill is that it should be able to reproduce the available 
constraints concerning the luminous sources. Whereas a relatively 
solid consensus has been reached on the luminosity function, spec- 
tra and evolutionary properties of intermediate redshift QSOs, some 
debate remains on the presence of yet undetected low-luminosity 
QSOs powered by intermediate mass black holes at high redshift 
(Ricotti & Ostriker 2004b). The evolution of stellar radiation is in- 
stead much less certain. Models corroborated by the observation of 
an excess in the Near InfraRed Background (NIRB) in the wave- 
length range 1^ fim (Salvaterra & Ferrara 2003; Magliocchetti, 
Salvaterra, & Ferrara 2003; Kashlinsky et al. 2004; Cooray et al. 
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2004) require that the first, metal-free (PopIII) stars were very mas- 
sive if they are to account for this otherwise unexplained excess. 
The same data require that a quite rapid transition from PopIII to 
"normal" PopII stars must have occurred at about a = 9, probably 
driven by an increase of the metallicity in the active cosmic star for- 
mation sites above the critical value Z^it = 1O-^±^Z0 (Schnei- 
der et al. 2002; Schneider et al. 2003). Even after that, predicting 
the joint reionization and star formation histories self-consistently 
is not an easy task as mechanical and radiative feedback processes 
can alter the hierarchical structure formation sequence of the un- 
derlying dark matter distribution as far as baryonic matter is con- 
cerned. 

Hence our present aim is to build up self-consistent models 
that can reconcile the discrepancy between the Gunn-Peterson op- 
tical depth and the first year WMAP observations, and at the same 
time satisfy a larger number of experimental requirements concern- 
ing the IGM temperature evolution. He reionization, the number of 
Lyman-limit systems, the NIRB excess interpretation, and the cos- 
mic star formation history. Although this might seem to a first sight 
a very ambitious aim, it turns out that once the relevant physics is 
included, these results allow a very clear scenario to emerge. 

Nothing comes for free, unfortunately, and there is a price to 
pay for this wealth of predictions that can be obtained via a rel- 
atively simple, although physically rigorous, semi-analytical ap- 
proach to this problem. When compared to the numerical simula- 
tions, we can treat many details of the reionization process only in 
an approximate manner (the shape of ionized region around sources 
and their overlapping, just to mention a few) and in terms of global 
averages (such as the filling factor and the clumping factor of ion- 
ized regions). However, we have to recall that in order to exploit the 
full power of the observational data available to constrain models, 
one must be able to connect widely differing spatial (and temporal) 
scales. In fact, it is necessary, at the same time, to resolve the IGM 
inhomogeneities (sub-kpc physical scales), follow the formation of 
the very first objects in which PopIII stars form (kpc), the radia- 
tive transfer of their ionizing radiation and the transport of their 
metals (tens of kpc), the build up of various backgrounds (Mpc), 
and the effects of QSOs (tens of Mpc). Thus, a proper modelling 
of the relevant physics on these scales, which would enable a di- 
rect and simultaneous comparison with all the available data set 
mentioned above, would require numerical simulations with a dy- 
namical range of more than 5 orders of magnitude, which is far cry 
from the reach of our current computational technology. To over- 
come the problem, simulations have typically concentrated on try- 
ing to explain one (or few, in the best cases) of the observational 
constraints. It is therefore difficult from these studies to understand 
the extent to which their conclusions do not conflict with a different 
set of experiments other than the one they are considering. 

The quest for a more general, data- supported scenario is the 
main motivation for this study. Although our semi-analytical ap- 
proach provides a way to achieve this aim in a very straightforward 
manner and, as we will see, with very satisfactory results, it is im- 
portant that its uncertainties are kept under control as much as pos- 
sible. For this reason we also compare its predictions with those 
from numerical simulations focusing on a more restricted aspect of 
the general picture. The final hope is that heavily data-constrained 
models, like to one presented here, will help us to identify the 
unique way in which cosmic reionization occurred. 

The paper is organized as follows: In the next few sections, we 
develop the analytical formalism for studying the ionization and 
thermal history of the IGM, implementing most of the essential 
physical processes. We discuss the evolution of the ionized regions 



in Section 2, taking into account the inhomogeneous density dis- 
tribution of the IGM. The explicit form of the IGM density distri- 
bution used in this work is discussed in Section 3. The estimation 
of the number of ionizing photons from different types of sources 
is given in Section 4. The formalism for studying the thermal and 
ionization histories of different phases of the IGM is discussed in 
Section 5. We discuss some properties of the absorption systems 
derived from our formaUsm in Section 6. Finally, we compare the 
model predictions with various observations and discuss the effect 
of varying different free parameters in Section 7. The final section 
presents the summary of our main results and compares the model 
with other semi-analytical models and numerical simulations. 



2 EVOLUTION OF IONIZED REGIONS 

In the following sections, we shall highlight the crucial points of 
our formaUsm and try to develop a consistent picture of cosmic 
reionization which derives from all known experimental constraints 
on such a process. Such formalism allows us to track (i) the evo- 
lution of the volume filling factors of ionized hydrogen and he- 
livmi regions and (ii) the thermal and ionization evolution in each 
of these regions separately and self-consistently. 



2.1 Volume filling factor of ionized regions 

Let us first start with individual ionized regions: these can either be 
singly-ionized hydrogen (HII) or doubly-ionized helium (Helll). 
We do not consider the analogous Hell regions, as typically these 
match the HII ones and, in any case, very sparse observational data 
have been collected concerning this ionization state of the IGM. 
To a first approximation, it is usually assumed that the IGM can be 
described as an imiform medium with small scale climipiness taken 
into account through a so-called "clumping factor". In such cases, 
the evolution of the volume filling factor Qhii for the HII regions 
will be described by (Shapiro & Giroux 1987) 

— t;^ — = VHiiCmi^afl (1) 

dr uh a-' 

where rtph is the rate of ionizing photons per unit volume, riH is 
the hydrogen atoms nvmiber density, Chii is the clumping factor 
for the ionized regions, a is the expansion factor, and an is the 
recombination coefficient. Under the assumption of homogeneity, 
the same equation holds for the ionized hydrogen mass fraction as 
well. One can write a similar equation for the volume filling fac- 
tor Qneiii of Helll regions too. In this picture the reionization is 
said to be complete when individual ionized regions overlap, i.e., 
Qi = 1. The above equation can be solved given a model for the 
source term hph, and a value for the clumping factor d. There are 
few points worth noting: (i) The quantity (corresponding to the re- 
combination term on the right hand side) QHiiCmirienHaR gives 
the average recombinations per unit time per unit volume in the 
universe, (ii) Also note that we have implicitly assumed that each 
photon is absorbed shortly after it is emitted, i.e., its mean free path 
is much smaller than the horizon size (which is true for z > 2). 
In case the photon is able to travel large distances before it is ab- 
sorbed, one has to take into account the fact that photons might 
be redshifted below the ionization edge of hydrogen, (iii) Finally, 
one should remember that the quantity hph takes into account only 
those photons which escape into the IGM. The number of photons 
produced by the source can be much larger; however, a large frac- 
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tion of those photons will be absorbed in ionizing the high density 
halo surrounding the luminous source. 

2.2 Inhomogeneous Reionization 

In the above picture, the inhomogeneities in the IGM are consid- 
ered simply in terms of the clumping factor in the effective recom- 
bination rate without taking into account the density distribution of 
the IGM. The importance of using a density distribution of the IGM 
Ues in the fact that regions of lower densities will be ionized first, 
and high-density regions will remain neutral for a longer time. The 
main reason for this is that the recombination rate is higher in high- 
density regions where dense gas becomes neutral very quickly. Of 
course, there will be a dependence on how far the high density re- 
gion is from an ionizing source. A dense region which is very close 
to an ionizing source will be ionized quite early compared to, say, 
a low-density region which is far away from luminous sources. 

The first effort in addressing such issues were carried out by 
Miralda-Escude, Haehnelt, & Rees (2000) (MHR, hereafter). To 
summarize the main point of this approach, consider, first, the sit- 
uation where all the individual ionized regions have overlapped 
(the so-called post-overlap stage; Gnedin 2000). In this scenario, 
all the low-density regions of the universe will be highly ionized, 
while there will be some high density peaks (like the collapsed sys- 
tems) which will still remain neutral. Thus, as a reasonable first 
approximation (MHR), we assume that all regions with overdensi- 
ties A < Ai are ionized (the index i refers to the different ionized 
species), while regions with A > Aj are completely neutral, with 
Aj increasing as time progresses (i.e., more and more high density 
regions are getting ionized). We shall discuss the equation govern- 
ing the evolution of Ai later in this section, while the results for 
different model parameters will be discussed in Section 7 (see, for 
example, Fig. 1). Note that, according to this scenario, the reioniza- 
tion is defined to be completed when all the regions with A < Aj 
are ionized - one does not need to ionize the whole IGM to com- 
plete the reionization process. The effect of this assumption is that 
only the low-density regions will contribute to the clumping fac- 
tor (regions whose density is above A, are assumed to be neutral, 
hence they need not be taken into account while calculating the 
clumping factor). 

The situation is slightly more complicated when the ionized 
regions are in the pre-overlap stage. At this stage, it is assumed that 
a volume fraction 1 — Qi of the universe is completely neutral (ir- 
respective of the density), while the remaining fraction of the 
volume is occupied by ionized regions. However, within this ion- 
ized volume, the high density regions (with A > Ai) will still be 
neutral. Once Qi becomes unity, all regions with A < Aj are ion- 
ized and the rest are neutral. The high-density neutral regions man- 
ifest themselves as the Lyman-limit systems in the QSO absorption 
spectra. The reionization process after this stage is characterized 
by increase in Ai - implying that higher density regions are getting 
ionized gradually. 

To develop the equations embedding the above physical 
picture, we need to know the probability distribution function 
P(A)dA for the overdensities. We shall discuss the form of ^'(A) 
in the next section, but given a P( A)dA, it is clear that only a mass 
fraction 

FM(Ai)= / "dAAP(A) (2) 
Jo 

needs be ionized, while the remaining high density regions will be 
completely neutral because of high recombination rates. The gen- 



eralization of equation (1), appropriate for this description is given 
by MHR (see also Wyithe & Loeb 2003) 

d[QHiiPM(AHii)] npi,{z) _ afl (T)nei?( Ahii) 
T7 = VHII ^ (3) 

where QmiaR{T)nmineR{Ami) gives the number of recombi- 
nations per unit time and volume. The factor i?(AHii) is the anal- 
ogous of the clumping factor, and is given by 

R{Ami) = dAA^ P(A) (4) 

Jo 

The reionization is complete when Qhii = 1; at this point a mass 
fraction Fm ( Ahii) is ionized, while the rest is (almost) completely 
neutral. 

Note that there are two unknowns Qhii and PM(AHn) in 
equation (3) [and similarly for the Helll regions]. For the post- 
overlap stage, we put Qi = 1, and can solve the equation for Ai. 
However, for the pre-overlap stage, we have to deal with both the 
unknown and it is thus impossible to solve it without more assump- 
tions. One can fix either Ai or Fm (the ionized mass fraction). 
There is no obvious way of dealing with this problem. In this work, 
we assume that Ai does not evolve with time in the pre-overlap 
stage, i.e., it is equal to a critical value Ac. To fix this Ac, we note 
that results do not vary considerably as Ac is varied from ~ 10 
to ~ 100. For definiteness, we take the value corresponding to the 
typical overdensity of collapsed haloes at the virial radius. The ex- 
act value of this overdensity depends on the density profile of the 
halo; it can be shown that it is ~ 59.2 for a isothermal profile, while 
it is « 63.7 for the NFW profile. In this paper, we can safely as- 
sume Ac = 60, and also assume that this critical value is the same 
for HII and Hein. Once Aj is fixed, one can follow the evolution of 
Qi until it becomes unity. Following that, we enter the post-overlap 
stage, where the situation is well-described by equation (1). 



3 PROBABILITY DISTRIBUTION P(A) 

There can be various approaches to determine the density distribu- 
tion of the IGM at various redshifts. It is clear that various compli- 
cated physical processes will not allow us to obtain a simple distri- 
bution from analytical calculations. One has to resort to either sim- 
ulations, or some sort of approximation schemes. In the approach 
followed by MHR, one uses a density distribution inspired by hy- 
drodynamical simulations. While such approaches are widely used, 
one should keep in mind that the limitations related to box size 
and resolution are inherent to every simulation. On the other hand, 
there are approaches where the density distribution is obtained from 
some approximation scheme. These approximation schemes seem 
to be reasonable in the Unear or quasi-linear regime of density fluc- 
tuations, while they are more inaccurate when non-linear effects 
creep in. 

In this work, we shall use such an approximation to describe 
the low-density IGM. There are several reasons to believe that the 
low-density regions of the IGM are well described by the lognor- 
mal distribution (see, for example, Choudhury, Padmanabhan, & 
Srianand 2001; Choudhury, Srianand, & Padmanabhan 2001; Viel 
et al. 2002), which is shown to be in excellent agreement with nu- 
merical simulations (Bi & Davidsen 1997). However, the lognor- 
mal distribution seems to fail at very high densities (say, when the 
densities are typical to that of collapsed haloes). To see this, note 
that at high densities (A S> 1), the lognormal distribution has the 
limiting form P(A) ~ ^^-inA/cr -i^ whereas it is expected that 
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it should follow a simple power-law -P(A) ~ A''. In fact, if dark 
matter haloes have a density profile of the form r'''^'''^^^\ then it 
can be shown that P(A) should fall as A''. For an isothermal pro- 
file, one has (3 = —2.5, while for a NFW profile, the value of /3 
varies from —4 in the center of the halo to —2 at the virial radius. 

We thus assume the probability distribution of the overdensi- 
ties to be lognormal at low densities, changing to a power law form 
at high densities: 



P(A) = ^l_e-ii^('"^-'')' 



if A < Av 



ctAvV^tt J 

(5) 

where a is rms linear mass fluctuations in baryons.^ If we want the 
derivative of P(A) to be continuous at the transition overdensity 
Av^, then it follows that the slope (3 and Ay must be related by 



1_ J_(lnAv'-M) =/? 



(6) 



The parameters A and ji are determined by normalizing the vol- 
ume and mass to unity. In most cases, the power law form of the 
probability distribution comes into effect only for the post-overlap 
phase. 

In general, haloes could have /3 varying with the distance from 
center depending on the density profile. Simulations also suggest a 
redshift evolution of /3, with (3 = —2.5 at high redshifts, while 
it is « —2.2 at redshifts around 2 (Miralda-Escude et al. 1996). 
In this work, however, we assume a constant value of —2.3 inde- 
pendent of the redshift and the distance from the center. This is 
reasonable as long as we are not too close to the center of the halo; 
also most of the results are insensitive to (3 as long it is within the 
range [—2.2, —2.5]). Given a value of P, the value of /S.v will, in 
general, evolve with time. 

To proceed further in the solution of equation (3), one has 
to estimate two quantities: (i) the photon production rate hph(z), 
and (ii) the temperature, T, of the ionized regions. We discuss the 
method adopted to obtain these estimates in the next sections. 



4 PHOTON PRODUCTION RATE 
4.1 Photons from galaxies 

One can use the Press-Schechter formalism to estimate the mass 
function of dark matter haloes of mass M which collapsed at a 

redshift Zc. In this paper, we shall use the formalism of Sasaki 
(1994) for calculating the formation and merging rates of dark mat- 
ter haloes. Assuming a model where the star formation rate (SFR) 
peaks around the dynamical time-scale of the halo, it will form stars 
at the rate (Chiu & Ostriker 2000; Choudhury & Srianand 2002) 



Msf(M,z,«c) = esF ( 



t{z)-t{Zc) 



(7) 



^ Throughout this paper we will assume a flat universe with total matter, 
vacuum, and baryonic densities in units of the critical density of Qrn = 



0.27, = 0.73, and ^5/*^ 



0.024, respectively, and a Hubble con- 



stant of Hq = 100 /i km s~i Mpc~i, with h = 0.72. The parameters 
defining the linear dark matter power spectrum are cs = 0.9, n = 0.99, 
dn/dlnA; = 0. 



where esF is the efficiency of star formation. We can then write the 
cosmic SFR per unit comoving volume at a redshift z. 



Psf{z) 



poo pc 

/ dzc 

J z J M 



dMMsF{M, z, Zc)N{M, z, Zc)m 



where N{M, z, Zc)dMdZc gives the number of haloes within 
a mass range (M, M + dM) formed within a redshift interval 
{zc,Zc + dzc) and surviving down to redshift z. The lower inte- 
gration limit, Mmin(2c), takes into account the fact that low mass 
haloes will not be able to cool and form stars. Also, note that we do 
not take into account possible disruption of star forming haloes due 
to energy injection from exploding supemovae. We shall discuss 
the choice of Minin(zc) in detail in a later section. 

Putting all the relevant expressions together, one can write the 
cosmic SFR in a neater form: 



Psf{z) 



where 



1 /"= 



dZcesF Fi{z, Zc)I{zc) 



(9) 



Fi{z,Za) 



and 



D{zc) 



D{zc)H{zc) 



D{zc) t{z)-t{z 



(1 + Zc) 



"-dyn 



{Zc 



(10) 



dv 



i'(M„i„[zcl) 



— e 



-1/^/2 



(11) 



In the above expressions, D{z) gives the linear growth factor of 
dark matter perturbations, and i/(M) — 5c/[D{z)a-DM{M)]. In 
this paper, we fix the critical overdensity for collapse, Sc, to 1.69. 

The rate of ionizing photons per unit voliune per unit fre- 
quency range would be 



n„,G{z) 



dN, 



dM 



PB 



D{z) 



dZcesFfescFl{z,Zc)T{Zc) (12) 



where /esc is the escape fraction of ionizing photons from the star 
forming haloes, and dA^,y/dM gives the number of photons emitted 
per frequency range per unit mass of stars. Given the spectra of 
stars of different masses in a galaxy, and their Initial Mass Function 
(IMF), this quantity can be computed in a straightforward way. The 
IMF and spectra will depend on the details of star formation and 
metallicity, and can be quite different for Pop II and Pop III stars. 
We shall return to this point at the end of this section. Given the 
above quantity, it is straightforward to calculate the total number 
density of ionizing photons emitted at a particular frequency range 
[t'min, t'max] by galaxics per unit time. 

4.2 Ionizing photons from PopII/PopIII stars 

While considering the number of photons from stars, one has to 
keep in mind that the nature of stars which were formed early 
could be very different from what we observe at lower redshifts. 
The physical motivation behind this assumption is the coimnonly 
accepted view that early stars are metal-free and hence the IMF 
would be top-heavy, i.e., dominated by high-mass stars (Schneider 
et al. 2002). Currently, the observational support for this assump- 
tion comes mainly from the analysis of the cosmic NIRB data (Sal- 
vaterra & Ferrara 2003; Magliocchetti, Salvaterra, & Ferrara 2003). 
Hence, in this work we consider the possibility that at early red- 
shifts, a population of mctal-free, massive (PopIII) stars produced 
a large number of photons. On the other hand, star formation at low 
redshifts, as we know, is dominated by high metallicity PopII stars 
with the usual Salpeter-Uke IMF. It is believed that the transition 
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from PopIII to PopII stars occurred because of an increase in tiie 
metallicity in thie active star-forming regions above a critical value 
ofZcrit = 10~^*^Zo. However, it is not clear how abrupt the tran- 
sition has been, due to the poorly known amplitude of the optical 
background at 0.6-0.8 fim. We shall therefore denote this transition 
redshift by «trans, and study the effects of its variation. We should 
mention that PopIII star formation may continue at lower redshifts 
with a decreasing rate, however the interpretation of the NIRB re- 
quires that their contribution to the ionizing flux be negligible with 
respect to the PopII stars. Although some metal pollution has to 
take place for transition to PopII stars, it is necessary that the ma- 
jority of the NIRB-contributing PopIII stars have masses such that 
a black hole is the final product of their evolution (Salvaterra & 
Ferrara 2003). 

From the above, it thus turns out that the star formation is 
made up of two components: 



Mph,G(z) = "ph, PopII (z) + nph,PopIIl(z); 



(13) 



the previous expression involves two free parameters which are the 
ionizing photon efficiencies of the Popn and Popni stars, namely: 
epopii = esF.Popii/osc, PopII and epopiii = esF.Popiii/esc.Popiii- 
We assvune that stellar population within the galaxies formed at 
z > Ztmns is dominated by PopIII stars, while at lower redshifts 
it is dominated by PopII stars. Since the SFR peaks at the dynam- 
ical time scale after the formation of the halo and decreases ex- 
ponentially thereafter, the formation rate of PopIII stars continues 
for some amount of time after Ztrans, and hence the transition is 
not instantaneous. By this assumption, we are implicitly neglect- 
ing the effects of galactic self-emichment on the IMF transition. 
Although locally (i.e., inside a given galaxy), the transition from 
PopIII to PopII star formation mode can occur at epochs different 
(either earlier or later) than ztrans, the studies of NIRB data set the 
epoch when the bulk of the cosmic stars became PopII at ^trans. 

To calculate the number of photons produced per unit mass of 
PopH star formed [dAv/(lM]popii, we use the stellar spectra cal- 
culated using STARBURST99 ^ (Leitherer et al. 1999), with metal- 
licity Z = 0.001 = 0.05Zg and standard Salpeter IMF. The num- 
ber of photons per unit mass of star formed, when integrated over 
appropriate frequencies, gives (8.05,2.62,0.01) x 10*'" Mq^ for 
(HII, Hen, Helll) respectively. For calculating [dA^„/dAf]popiii, 
we first assume that the IMF of the PopIII stars is dominated by 
very high mass stars. Since the number of photons produced per 
imit mass of stars formed is somewhat independent of the stellar 
mass for high mass stars, it is independent of the precise shape 
of the IMF (as long as the IMF is dominated by high mass stars). 
The value of [dA^^/dM]popiii is calculated using stellar spectra for 
high mass stars (> SOOM©) as given by Schaerer (2002). When 
integrated over appropriate frequencies, this quantity is equal to 
(2.69, 4.46, 1.36) x lO^^M"^ for (HII, Hell, HeHI) respectively. 

4.3 The minimum mass for the star-forming haloes 

The quantity Afmin(zc) in equation (8) depends on the cooling ef- 
ficiency of haloes. We assume that molecular cooling is active at 
high redshifts, and M^in{zc) increases with decreasing z (Fuller 
& Couchman 2000). However, it is also possible that the molecu- 
lar cooling within the mini-haloes could be highly suppressed due 
to photo-dissociation of hydrogen molecules. We shall study the ef- 
fect of mini-halo suppression in Section 7.5. At lower redshifts (say 



z < 10), when the mass function of collapsed haloes is dominated 
by relatively high-mass haloes, the SFR is more or less independent 
of whether molecular cooUng is effective or not. 

There is a further factor which needs to be taken into account 
- the radiative feedback from stars. Once the first galaxies form 
stars, their radiation will ionize and heat the surrounding medium, 
increasing the mass scale (often referred to as the filtering mass) 
above which baryons can collapse into haloes within those regions. 
The minimum mass of haloes which are able to cool is thus much 
higher in ionized regions than in neutral ones. Since we are consid- 
ering a multi-phase IGM, one needs to take into account the heating 
of the ionized regions from the beginning (even before the actual 
overlap has started). As we compute the temperature of the ionized 
region self-consistently, we can calculate the minimum circular ve- 
locity of haloes that are able to cool using the relation: 

2A;boitzT 



2 



(14) 



where /i is the mean molecular weight: /imp = pb/nt- Typically, 
for a temperature of 3 x W^K, the minimum circular velocity is ~ 
30-50 km s~^, which is comparable to the filtering scale obtained 
from numerical simulations of Gnedin (2000) after the universe has 
reionized. 

Note that the above value of Vc will evolve according to the 
temperature of the gas. An alternative feedback prescription is to 
fix the minimum circular velocity of haloes allowed to form stars 
at a given value, which can be taken to be in the range Vc = 35- 
50 km s"^ (Gnedin 2000; Kitayama et al. 2001). We shafi explore 
the effect of such different feedback prescriptions later in Section 
7. 



4.4 Photons from QSOs 

In order to calculate the nimiber of ionizing photons from QSOs, 
we shall follow the simple formaUsms developed by Wyithe & 
Loeb (2003) and Mahmood, Devriendt, & Silk (2004). The only 
difference in our approach is that we use a different prescription for 
calculating the merging and formation of dark matter haloes. In the 
previous works, the merging rates of dark matter haloes were based 
on the formalism by Lacey & Cole (1993), while our model uses 
the Sasaki (1994) methodology. For the reason of completeness, 
we include the details of the model for calculating the limiinosity 
function of QSOs in this section. 

The key assumption to calculate the luminosity of QSOs is 
that the mass of the accreting black hole Mbh is correlated with the 
circular velocity Vc of the collapsed halo through the relation: 



Mbh = ev° 



(15) 



It is argued that a = 5 in a self-regulated growth of super-massive 
black holes (Silk & Rees 1998), a value found to match observa- 
tions of local universe. 

It is then reasonable to assume that the black hole radiates with 
the Eddington luminosity (in the B-band), given by 



Lq,b 
which gives 



5.7 X 10' 



0,s 



Mn 



a/3 



(16) 



(17) 



http://www.stsci.edu/science/starburst99 



where we have used the relation between the circular velocity and 
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the mass of a virialised halo (Choudhury & Padmanabhan 2002) to 
obtain 



5.7 X lO^olO"*" 



H2(2)Avir(^) 



i/6 



1 - 



3i?2(z)Avir(^). 



a/3 



with 



£0 = 



£(159.4 km s-i)° 



(18) 



(19) 



Note that, even if the black hole radiates in a sub-Eddington rate, 
the corresponding uncertainty can be absorbed into the value of /?. 
The luminosity function of QSOs (i.e., the number of QSOs per 
unit comoving volume per unit luminosity range) will be given by 



ip{LB,z)dLB 



f 



dZcN{M,z,Zc)dM 



(20) 



where M and Lb are related by equation (17). If we assume that 
each QSO lives for a time tqso(2) ^ (d/a)"^, then the QSO activ- 
ity can be taken to be almost instantaneous and we can approximate 



dz 



I dzc N{M, z, z, 
Then 

^P{Lb,z) = '^t^,^{z)N(M,z,z) 



—t^,o{z)N{M,z,z) 



dM 



(21) 



3 



MNm{z)v 



dlB 

D{z) 



D{z)H{z) 



H{z)tqBo(z) 

(22) 



We fix tqso(2) = 0.035tdyn(z) (Mahmood, Devriendt, & Silk 
2004). Given a = 5, one can fix the free parameter eo by com- 
paring the model with observations of QSO luminosity function. 

This simple procedure works very well at intermediate and 
high redshifts, but fails to match the observations at low redshifts. 
One can introduce a phenomenological function in equation (22) to 
take into account the break in the luminosity function at high lumi- 
nosities. We find that a modified iI>{Lb,z) of the form (Mahmood, 
Devriendt, & Silk 2004) 



ijj{LB,z) — > iP{Lb,z) exp 



M 



1011.25+^M( 



© 



(23) 



is good enough in matching the low redshift observations. This cut- 
off is at very high luminosities and has virtually no effect at 2; > 3. 

Given the luminosity function, the rate of ionizing photons 
from QSOs per unit volume per unit frequency range will be 



(«) = / dLB-ip{LB,z) 



L^Lb) 



(24) 



We next use the mean UV QSO spectrum (Schirber & Bullock 
2003) 



L.{Lb) 



-'B j^qIS.Oo 



ergs s~iHz i'O.s 
which then gives, 

^10i«°5ergss-iHz 



(25) 



f 

Jo 



Lq,b 
dl/s Lb iP{Lb,z) 



hpVH 



(26) 



The rate of ionizing photons from QSOs is obtained simply by in- 
tegrating over all relevant frequencies. 

This simple phenomenological model suits very nicely the 
purposes of this paper. Alternatively, one can simply use the ob- 
served luminosity function for QSOs for calculating the number of 
ionizing photons. Since the above model matches the perfectly well 
with observations at low redshifts, none of the results would be af- 
fected. At high redshifts, say z > the contribution from QSOs 
are negUgible compared to that of galaxies, and can be ignored. 



5 THERMAL EVOLUTION 

In the previous sections, we have discussed the evolution of the 
ionized volume filling factors, and various physical quantities re- 
lated to them. It is clear that the baryonic universe will behave as 
three-phase medivun constituted by: (i) completely neutral regions, 
(ii) regions where hydrogen is ionized and helium is singly ionized, 
and (iii) a region where both species are fully ionized. We have as- 
sumed that the ionization front for the doubly-ionized helium can 
never overtake that for ionized hydrogen - which is found to be 
always true for the type of spectra we are using for the ionizing 
sources (note that a much harder spectrum can always make the 
ionization front for the doubly-ionized helium leading the ionized 
hydrogen region). 

The thermal evolution equations are solved separately for each 
of the three regions. In the absence of heating sources the evolution 
is nearly trivial for the neutral region, with the temperature decreas- 
ing adiabatically. However, there is always a background of hard 
photons which can heat the neutral IGM. Since the temperature of 
the neutral region hardly affects the reionization history, we ignore 
such hard photons and let the temperature of the neutral regions de- 
crease adiabatically. In the ionized regions, the temperature can be 
calculated using the dynamical equation 



where 

Xi 



dE 



dt 3feboitznB(l -I- ^;)^ dt 



(27) 



PB 



(28) 



and dE/dt gives the net heating rate per baryon. For most pur- 
poses, it is sufficient to take into account the photoionization heat- 
ing and recombination cooling (and Compton cooling off CMB 
photons, which can be important at high redshifts). For example, 
in the regions where only hydrogen is ionized, we have 

'^^ = nm{l + z)^ f dv \h{z; v) j^"^^)^^ a_g(i/)/ip(z/ - vh) 



dt 



Qmi{z) 



- R{Ami)aRciT)minneil + zf 



(29) 



where \h(z\ u) is the proper mean free path for hydrogen ion- 
izing photons with frequency v > vh- It is found from obser- 
vations at low redshifts that Xh{z; v) oc f^'^; this is understood 
from the frequency-dependence of the absorption cross section 
(Jh{v) oc v^'^ , and the column density distribution of QSO ab- 
sorption Unes dN/dNm oc N^^^^^ (Petitjean et al. 1993). We 
assume this relation to be valid for all redshifts. Note that this 
frequency dependence of the mean free path hardens the ioniz- 
ing spectrum. However, at frequencies below the ionization thresh- 
old of Hell, the diffuse recombination radiation from the IGM 
tends to compensate for this hardening (Haardt & Madau 1996). 
Given this, we assume \h{z;ii) = \Hfi{z) for v < z/Heii and 
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\h{z;v) = \H,oiz){i^/meiiy'^ for u > meii- The procedure 
for calculating Xnfiiz) = Xh{z; v = vh) is described in the next 
section. 

The equation for evolution of the temperature has to be sup- 
plemented by those for the ionization of the individual species. In 
the most general case, one has three independent species Xi = 
{Xhi, Xnei, ^Heiii}, with other species being given by 

Xhii = 1 — Y — Xhi; 
Y 

XHell = -J — Xnel — XHelli; 

Xe = Xhii + Xhoii + 2XHeiii (30) 

where Y = 0.24 is the helium weight fraction. For example, the 
evolution of Xhi in the HIl region is given by 

dXHi 



dt 



-Xhi(1 + z 



(^mi(z) 



+ R{Ami)aR{T)XmiXe — {l + zf (31) 

nip ' 

Similar equations, though slightly more complicated, can be writ- 
ten down for other regions too. 

In passing, note that the volume emissivity is given by 



ev(z) — hv{z)hpiy{l + z)'^ 



(32) 



while the ionizing flux for a particular species i [one of HI, Hel or 
Hell] is given by 



-n„(z)hpu{l + z) 



4n 47r 
The photoheating rate for a particular species is given by 

J. 



(33) 



r(2) = 47r 



f 



dv^^ai{v)hp(y - i^min) 
hpv 



du Xi{z; v)nu{z)ai{v)hp{v — z/„ 



(34) 



where Vmin is the threshold frequency for the species considered. 
The photoionization rate is given by 

J. 



Tpiiz) = 47r 



di^-^ai{v) 
hpv 



= (1 + 4 



f 



di/ \i{z; v)h^[z)ai{v) 



(35) 



Note that we have included the clumping term J?(Ai) in the 
expressions (29) and (31). As more and more regions of higher 
densities get ionized, the value of ii(Ai) becomes larger which, 
in turn, gives a larger value of the temperature. Thus the temper- 
ature T of the ionized regions obtained from the above system of 
equations are essentially weighted by the mass of the correspond- 
ing regions. In this sense, one can assume T to be an estimate of 
the mass-averaged temperature of the region. We should emphasize 
that T is not the rigorously-defined mass-averaged temperature, but 
should be considered as a simple approximation in the ionized re- 
gions. If the quantities T and Xi defined in this section are approx- 
imate estimates of the mass-averaged values in the ionized region, 
then the global mean values Tgiobai and Xgiobai.i can be obtained 
by weighted averages over different regions, according to the mass 
fraction of the corresponding region. Also, note that the above T is 
not same as the conventional To, which is defined as the tempera- 
ture of gas at the mean IGM density (A = 1). Similarly, the ion- 
ization fractions defined above need not correspond to the values at 
the mean density. The values at the mean density (i.e.. To, Xhi.o) 



can be solved using the same equations (27) and (31), but without 
putting in the clumping factor i?(AHii), i.e., 

— --2H(z)To-^-^^^ 

nm.OT^ — 7-T - aHc(-ro)nHii,one,o(l + z) 



QhiiC^:) 



(36) 



and 

dX^ifi 



= -Xhi.o ^^'^,^\ +Qi{(ro)XHii,oXe,o^(l-h-?)^37) 
dt (3hii(z) rup^ ' 



5.1 Mean free path for photons 

The mean free path for ionizing photons depends on the size and 
topology of the ionized regions. Hence, for a self-consistent cal- 
culation of the mean free path, one has to use the evolution of the 
volume filling factor of the ionized regions. Note that we only have 
statistical information about the fraction of volume and mass which 
is ionized, i.e. we do not calculate the size of individual ionized re- 
gions. 

Given this situation, we use a simple model developed by 
MHR to calculate the mean free path Xi^o(z) for photons (at 
v = Vh). As discussed in MHR, their method is a good approx- 
imation when a very high fraction of volume is ionized. It is clear 
that a photon will be able to travel through the low density ionized 
volume 



Fv{Ai) 



Jo 



dA P(A) 



(38) 



before being absorbed. In the simple model, one assumes that the 
fraction of volume filled up by the high density regions is 1 — Fy, 
hence their size is proportional to (1 — FvY^^, and the separation 
between them along a random line of sight will be proportional to 
(1 — Fv)~^^^, which, in turn, will determine the mean free path. 
Then one has 



Xi,o{z) = 



Ao 



[l-Fv(A0]2/3 



(39) 



where we can fix Ao by comparing with low redshift observations. 
In fact, it has been suggested (from simulations and structure for- 
mation arguments; MHR) that Ao should be determined by the 
Jeans length which, in turn, is determined by the minimum circular 
velocity for star-forming haloes [equation (14)]: 



Xb{z) = Ho 



7 



3n,n(i + z)' 



(40) 



where 7 is the adiabatic index. In this work, we assume Ao oc 
Xb(z), with the proportionality constant being determined by com- 
paring with low redshift observations. The mean free path for pho- 
tons at 1/ = Vh is given by the typical separation between the 
Lyman-limit systems, which is observed to be ~ 33 Mpc at 2; = 3. 
From the knowledge of Ai,o(^) on can then predict the number of 
Lyman-limit systems per unit redshift range through the relation 
(Madau, Haardt, & Rees 1999; Miralda-Escude 2003) 



dA^LL _ c 

dz ~ y/^\ifi{z)H{z){l + z) 



(41) 



which can be directly compared with available observations at 2 < 

Z <4:. 
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6 PROPERTIES OF ABSORPTION SYSTEMS 

In this section, we discuss about how to obtain some of the prop- 
erties of the IGM when they are observed in the absorption spectra 
of QSOs. 

Given the probability distribution, the estimates of the tem- 
perature, T, and of the neutral hydrogen density nui, we can com- 
pute the mean transmitted flux as will be observed in the absorption 
spectra of QSOs. The optical depth at a given point is given by 

r(x,^;) =/„nHi(x,2)(l-F0f (42) 

where la = 4.45 x 10^^*^ cm^"^ is a constant. It is natural to 
assume that the matter and radiation in the IGM are in photoion- 
ization equilibriimi; in that case, the neutral hydrogen density is 
related to the baryonic overdensity through 

nHi(x,«) = nHi,o(^)A'-^-°''^(x,2) (43) 

The optical depth in the ionized region will then be 

t{x,z) = ^(z)A^-^-°''^(x,«) (44) 

where 

A{z) = 7anHi.o(^)(l + ^fj^^ (45) 

In general, one should use the global average value of nHi,o(^) in 
the above expression, taking into account the neutral and different 
ionized regions. The transmitted flux is 

F(x,«) = e-^(-'^); (46) 
with its global mean given by 

F{z) = dA e-^(")^ P(A) (47) 

Jo 

The equation of state (EOS) can be written in terms of the adiabatic 
index as T oc A''^^; note that sometimes, and somewhat confus- 
ing, 7 rather then 7 — 1 is defined as the slope of the temperature 
- density relation. It is, in principle, possible to compute the value 
of 7 by studying the evolution of the temperature for fluid elements 
of different densities. However, such a study is somewhat beyond 
the scope of this paper, and will be reported somewhere else. As far 
as this work is concerned, we keep 7 as a free parameter varying 
in the range 1 < 7 < 2.4 and compute the transmitted flux over a 
wide range of values of 7. Note that the choice of 7 does not affect 
any of our other results. 

Another quantity that can be readily estimated from our mod- 
els is the optical depth of CMB photons due to Thomson scattering 
with free electrons. This can be written as 

,z[t] 

Te\{z) = arc / dt nglobal.e (1 + z) (48) 

Jo 

where rigiobai.e is the global average value of the comoving elec- 
tron density. We neglect additional small contributions to Tei aris- 
ing from early X-ray sources but we do include reUc free electrons 
from cosmic recombination (Venkatesan, Giroux, & Shull 2001). 



7 RESULTS 

In this section, we present the main results for our model along 
with their interpretation. In the first part, we analyze what we con- 
sider the "fiducial" model in terms of the choice of free parameters, 
and confront it with observations. We shall see that this "fiducial" 



model seems to match all the available experimental data, thus jus- 
tifying our analysis. 

The two main free parameters of the model are the ionizing 
efficiencies of the two stellar populations, both of which are quite 
uncertain. For "normal" Popll stars, the fiducial values are taken 
to be esF,Popii = 0.1, which is partly constrained by low red- 
shift observations of the cosmic SFR (Nagamine et al. 2004), and 
/esc = 0.5%, yielding epopii = 0.05%. The estimates of the es- 
cape fraction are quite uncertain, and we shall discuss the effects of 
its variation later on. 

For the metal-free PopIII stars, the best observational con- 
straints come from the excess in the cosmic NIRB data (Salvaterra 
& Ferrara 2003; Magliocchetti, Salvaterra, & Ferrara 2003), which 
imply an upper limit on the combination esF,Popiii(l — /esc) = 
esF, PopIII — epopiii < 0.1 for top-heavy IMF. However, this is not 
sufficient to constrain epopiii. It is found that for low mass haloes 
(M ~ 10® M0), the escape fraction can be as high as /esc ~ 0.95 
(Whalen, Abel, & Norman 2004), while it can be significantly low 
for high-mass haloes (Ricotti & Shull 2000; Kitayama et al. 2004). 
Since the mass function of collapsed haloes is dominated by low- 
mass haloes at high redshifts, it is reasonable to use high values 
of /esc, particularly before reionization. Once the IGM is ionized, 
the star formation in low mass haloes is suppressed and hence one 
should probably use a lower value for the escape fraction appropri- 
ate for haloes of higher mass. Similarly there are no strong (theoret- 
ical or observational) constraints on the value of esp.popiii- Hence, 
for simplicity, we ignore any redshift or mass-dependence of /esc 
and esF, PopIII and use the fiducial values esp.Popiii = 0.8% and 
/esc = 90% so that epopiii = 0.7%. Note that all our results (ex- 
cept the SFR) are sensitive only to the combination esp.popiii/esc. 
Hence one can vary the individual values of the two parameters 
keeping the value of epopiii same, and still obtain identical results 
for the reionization history. For example, low-mass haloes might 
have /esc « 1; this impUes that our fiducial model requires a star- 
forming efficiency of esp.Popiii = 0.7%. On the other hand, in the 
ionized regions (where the mass function is dominated by higher 
mass haloes), the escape fraction could be significantly lower (say, 
/esc 10%); consequently, the implied value of the star-forming 
efficiency in our fiducial model would be esp.Popiii = 7%, which 
is probably higher than the theoretical expectations, but still does 
not violate the upper limit from NIRB data. However, if the escape 
fraction is significantly lower than this, it can be quite difficult to 
match the NIRB data and reionization constraints simultaneously. 
The effects of varying epopiii will be discussed later. 

Finally, the value of the transition redshift ztmns is taken to be 
10, which seems to be favoured by NIRB observations. The best- 
fit value 2;trans = 8.8 (Salvaterra & Ferrara 2003) obtained for a 
burst-Uke mode of SFR might be somewhat larger when we allow 
the PopIII SFR to decrease exponentially even after z = ztrans. 
As explained, radiative feedback, setting the minimum mass of the 
star forming haloes in the ionized regions, is implemented in our 
formalism. We shall discuss its effects in detail. 

This fiducial model, which seems to fit most of the observa- 
tions, is shown in Fig. 1. It also allows us to build a self-consistent 
reionization scenario, whose different predictions are compared to 
the available data in each of the 12 Panels composing the figure. 

The reionization history is well exemplified by the evolution 
of the filling factor of ionized regions Qi [Panel (b)]. According to 
the curves shown, hydrogen reionization must have taken place at 
redshift z fa 15 while the Hell reionization is completed around 
z w 12. We stress here that when Qi becomes unity, the regions 
having densities less than A, are completely ionized thus signi- 
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l+z 1+z 1+z 

Figure 1. The fiducial model which matches all available observations. Adopted parameters are epopiii = 7 X 10""^, epopii = 5 X 10~*, ztrans = 10. 
The Panels show as a function of redshift: (a) critical overdensity for reionization, (b) filling factor of ionized regions, (c) effective clumping factor, (d) specific 
number of Lyman-hmit systems, (e) ionizing photons mean free path, (f) mass-weighted temperature, (g) mean-density gas temperature, (h) cosmic star 
formation history, (i) photoionization rates for hydrogen, (j) photoionization rates for helium, (k) Gunn-Peterson optical depth, (1) electron scattering optical 
depth. See text for detailed description of different Panels. 



Table 1. Parameter values used in different figures throughout the paper. 



Fig. 


eSF.PopII 


/csCjPopII 


eSF.PopIII 


/csc.PopIII 


2^trans 




1 


0.1 


0.005 


0.007 


0.9 


10.0 


using equation (14) 


2 


0.1 


0.0005 


0.007 


0.9 


10.0 


using equation (14) 


3 


0.1 


0.1 


0.007 


0.9 


10.0 


using equation (14) 


4 


0.1 


0.005 


0.002 


0.9 


10.0 


using equation (14) 


5 


0.1 


0.005 


0.007 


0.9 


11.4 


using equation (14) 


6 


0.1 


0.005 


0.007 


0.9 


10.0 


50kms-i 
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fying the completion of reionization, while regions with higher 
densities are completely neutral. Also note that Qhii > QhoIii, 
thus showing that the Helll ionization front never overtakes the HII 
front. Interestingly, the evolution of QHeiii is markedly affected 
by the hydrogen reionization; there is a decrease in Qiidii when 
Qhii becomes unity. The reason for this is that as more and more 
regions get ionized, the feedback effect depresses Popin star for- 
mation, thus decreasing ionization rates. However, while this de- 
crease is marginal with respect to the wealth of hydrogen ionizing 
photons, it profoundly affects the emissivity of photons above 4 
Ryd. A more remarkable difference in the evolution of the filling 
factor QHeiii is however seen at redshifts 3.5 < ^ < 7. During 
this epoch the model predicts Helll recombination followed by a 
second Hell reionization occurring at « « 3.5; QHeiii drops to 
a minimum value of 0.4, i.e. He is primarily in the singly ionized 
state at « w 5 — 6. 

This behavior can be understood from Panels (i)-(j), where the 
ionizing rates of the sources responsible for H and He ionization are 
shown. From these Panels we see that at high redshifts (> 10), the 
ionizing flux is totally dominated by PopIII stars, which however 
fade off at lower redshifts. PopIII stars have a hard spectrum, as 
seen by comparing the values of H and He ionizing rates at high 
redshifts in Panel (i) and (j) respectively. It then follows that the 
PopIII stars can ionize Hell quite efficiently. Once their formation 
is quenched, there is little source of Hell ionizing photons until 
the QSOs take over the production of photons above 4 Ryd around 
redshifts of 5 [see Panel (j)]. Hence while the first ionization of 
Hell is controlled by PopIII stars, the second one is induced by 
QSO radiation. The transition epoch has been fixed in the fiducial 
model to ztrans = 10, and its effect on the rates can be clearly 
appreciated; also note that transition from PopIII to PopII stars is 
not instantaneous. Similar conclusions can be drawn by inspecting 
the star formation history in Panel (h): PopIII stars produce a first 
maximum in the star forming activity at z « 15 where psF « 
0.003Mo yr~^ Mpc~^, followed by a drop because of radiative 
feedback from ionized regions and a subsequent raise due to the 
increasing contribution of PopII stars leading to a less pronounced 
peak, PSF ~ Q.IMq yr^^ Mpc^'^ at z = 4. According to these 
results, only 0.3% of the stars produced by z = 2 need to be PopIII 
stars in order to achieve the first reionization. 

Reionization proceeds from regions of low density to over- 
dense ones, as shown by the evolution of the critical density Aj 
for the HII and Helll regions in Panel (a) . We recall that regions 
having densities above Aj are neutral, while a fraction Qi of the 
regions with densities lower than Ai are ionized. By 2 = 5 ioniza- 
tion fronts have been able to penetrate inside quite dense regions, 
with A > 10^, leaving therefore tiny islands of neutral gas (mostly 
in the vicinity or part of galaxies) in a sea of ionized plasma. As 
the HII ionization front always leads the Helll front, it is obvious 
that the value of Ahii is always higher than AhbIii. Note that Ai 
is constant as long as Qi < 1. Panel (c), illustrating the evolu- 
tion of the clumping factor i?(Ai) of the ionized regions, points 
towards the same behavior. Its evolution is essentially determined 
by the lognormal distribution, being essentially constant in the pre- 
overlap stage where Ai is constant. The corresponding evolution 
of the mean free path \i{z) for H- and Hell-ionizing photons is 
shown in Panel (e). 

The derived scenario can now be directly compared with a set 
of different experimental constraints. We have already mentioned 
that the value of esp.Popii is chosen such that the predicted SFR 
agrees with observations. However, it is still interesting to note 
from Panel (h) that our model reproduces the observed evolution 



of the SFR. We next consider the number of Lyman-limit systems 
per unit redshift range [Panel (d)], as computed from the mean free 
path for the hydrogen ionizing photons. The data points with error- 
bars are taken from Storrie-Lombardi et al. (1994). Another useful 
comparison involves the mean Lya transmitted flux, or, the effec- 
tive (Gunn-Peterson) optical depth [Panel (k)], defined as — In F, 
where F is the mean transmitted flux. The three curves, from right 
to left are for 7 = 1.0, 1.7, 2.4, respectively. Although the general 
raising trend shown by the data (Songaila 2004) is reproduced quite 
well by the model, there is some indication that 7 may vary, not un- 
expectedly, with redshift. In fact, a better agreement with the data 
is found if the adiabatic index increases from 7 = 1 (isothermal 
EOS) at z ~ 3 to 2.4 at z = 6. This behavior is opposite to simple 
expectations based on the assumption that reionization takes place 
at z ~ 6. Although is it true that at reionization the gas tends to ap- 
proach isothermality, afterwards evolving to higher more adiabatic 
state, we have to keep in mind that our fiducial model predicts that 
hydrogen reionization has occurred long before z = Q and there- 
fore it is physically plausible that at later evolutionary stages the 
EOS deviated from 7 = 1; this value is eventually recovered when 
Hell reionization occurs at z = 3.5. 

As an additional test of the model we next consider the IGM 
temperature. In particular. Panel (f) shows the estimates of the 
mass-averaged temperature for HII and Helll regions and the global 
one (we recall that this quantity is obtained by weighted averages 
over different regions, according to the mass fraction of the cor- 
responding region). Note that this quantity continues to increase 
with decreasing z at lower redshifts (post-overlap stage) due to the 
fact that the mass-averaged temperature is dominated by regions 
of higher densities, and as such higher density regions get ionized 
gradually, the mass-averaged temperature rises. As already men- 
tioned, the quantity usually derived from QSO absorption line ex- 
periments is the temperature at the mean density. To [Panel (g)]. 
Since To represents the thermodynamic state of the medium at 
mean density, its behaviour is markedly different from that of the 
mass-averaged temperature, particularly in the post-overlap era. As 
we mentioned before, the mass-averaged temperature rises because 
regions of higher densities get ionized, while the ionization of such 
regions do not affect To at all. Following the first reionization. To is 
boosted to 1.5 X 10* K; from there it decreases because of adiabatic 
expansion; however, the overall trend is relatively flat. A more pro- 
nounced decrease is prevented by the newly available Hell atoms 
from recombination of this species, which provides extra photo- 
heating to the gas. The curve seems to fit reasonably well the data 
(taken from from Schaye et al. 2000), although a few of them are 
more than l-cr away. Apparently, our model would better reproduce 
the experimental trend if the second Hell reionization could be de- 
layed by roughly 0.5 redshift units. However, note that the data is 
consistent with the fact that He is singly ionized (HII regions) at 
z > 3.5, while it is consistent with He being doubly ionized (Helll 
regions) at 2 < 3.5. Unlike the data points, the global tempera- 
ture in our models rises gradually rather than showing a sudden 
jump. This is related to the fact that QHeiii has a rather gradual 
rise. It might be possible that the data show a sudden jump in To 
because it is based on a few number of lines of sight. Since we ex- 
pect that there will be large fluctuations in the temperature along 
various lines of sight when Qi < 1, one thus has to compute the 
temperature along numerous lines of sight to get the global mean, 
which can then be compared with the model. On the other hand, 
our model can be used to generate different lines of sight with dif- 
ferent conditions, and can be compared with existing observations. 
In this regards, we should mention that there are other analyses of 
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the same data (Ricotti, Gnedin, & Shull 2000) which do not show 
the above mentioned sudden jump in To around z ~ 3.5 mainly be- 
cause of larger error-bars. In this sense, one requires much robust 
analysis of the observational data for more concrete conclusions. 
Finally, we turn in Panel (1) to the electron scattering optical depth, 
Tei. The gray shaded region is the l-cr limit as obtained from the 
first year WMAP data. The predicted values match quite well with 
the WMAP constraint. 

We conclude that our fiducial model can explain in a self- 
consistent manner and simultaneously all the available data exist- 
ing on cosmic reionization, showing also a strong predictive power, 
in spite of its simple and somewhat idealized implementation. The 
emerging picture is one in which the universe has been initially 
reionized by massive PopIII stars both in H and He; the subsequent 
disappearance of such exotic stars, required by the NIRB, caused 
Helll recombination, followed by its second reionization induced 



by QSOs. This evolution has produced little effect on hydrogen, 
which remained essentially ionized throughout, as its ionization 
state was maintained by normal PopII stars. 

In the following sections, we investigate different flavors of 
such pictures produced by variation of the free unknown model pa- 
rameters. 

7.1 Constraints on the PopII ionizing efficiency 

The efficiency parameter epopii for the low-redshift, normal stars 
is one of the most ill-determined parameters of our model. We reit- 
erate that we have, mainly for simplicity, assumed the parameter to 
be independent of z and galaxy mass, although there are hints from 
the analysis of the Sloan Digital Sky Sun'ey data that star forma- 
tion efficiency might vary with the halo mass as M^^^ for galactic 
stellar masses < 3 x 10^° M© (Kauffmann et al. 2004). 
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As discussed earlier, epopii is a the product of the star-forming 
efficiency esF.Popii and the escape fraction /csc.Popii of photons 
from the halo, both of which may vary with redshift and galactic 
mass. The first quantity esp.Popii can be reasonably constrained by 
comparing the model with observations of the SFR at lower red- 
shifts. In most cases, it is found to be ~ 10%. The largest uncer- 
tainty comes from the escape fraction /esc,Popii. It is quite difficult 
to determine it observationally, and the current limits from different 
kinds of modelling vary from a few per cent to about 50%. 

In the fiducial model explored so far, we adopted as an edu- 
cated guess epopii = 0.05% while Fig. 2 was devised with the aim 
of determining a lower limit on epopii from our models. However, 
one can realize from there that such determination is quite difficult 
to make from low {z < 6) redshift observations. The reason is that 
once epopii is decreased (in this case by 10 times with respect to 
the fiducial one) PopII stars play a sub-dominant role and observa- 



tions can be explained just with the ionizing photons from PopIII 
stars and QSOs. This does not necessarily mean that the contribu- 
tion from the PopII stars is insignificant or unnecessary, but it is just 
that the present observations have very little dependence on PopII 
star efficiency once below a certain threshold. We believe lower 
limits on epopii should be derived from some other considerations. 

On the other hand, it is possible to obtain a stringent upper 
limit on epopii, as shown by Fig. 3. In fact, if epopii > 0.01, then 
the number of ionizing photons produced will be too high, yielding 
a too low Lya optical depth (or equivalently, a higher transmitted 
flux) than the observed value at redshifts around ~ 3. Note that for 
a star forming efficiency of 10%, this upper limit corresponds to an 
escape fraction of 10%, which should be considered as quite strin- 
gent. Also to be noted is that the temperature of the Helll regions is 
lower than that observed around z ~ 3.5, hence this higher value of 
epopii is unable to match the high Tq at those redshifts. The reason 
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Figure 4. Same as in Fig. 1, but for epopiii = 0.002. 



for this is that with such a high production of photons from PopII 
stars, the ionizing spectrum is no more dominated by QSOs and 
thus effectively becomes much softer. This, in turn, results in lower 
values of To for Helll regions. 

There is one more point which needs to be clarified in this 
section. We have noted that the change in the value of epopii af- 
fects the low redshift physics only. On the other hand, one can see 
by comparing Panel (i) and (j) of Figs 1, 2 and 3 that the values 
of Fpi are lower at high redshifts when we increase epopii. This 
apparent contradiction arises from the fact the value of the photon 
mean free path Ao in equation (39) is chosen so as to match the 
low redshift observations. When epopii is increased, the number of 
photons is higher, leading to potentially larger ionized regions. In 
order to control the size of the regions so that they match low red- 
shift observations, we have to decrease the value of Ao. In turn, this 
decreases the ionization (and heating) rates at higher redshifts. 



7.2 Constraints on the PopIII ionizing efficiency 

As in the previous section, let us now try to determine the limits 
on epopiii (without changing the value of ztrans)- In contrast to the 
case for epopii, the major uncertainty in epopiii comes from the 
star forming efficiency esp,popiii- The escape fraction at such high 
redshifts is usually quite high (nearly 100%) for low mass haloes 
because typically only a few PopIII stars are required to completely 
ionize the parent galaxy (the masses of of PopIII-hosting galaxies 
are low). However, our understanding of the star formation process 
and constraints of parameters at high redshifts are quite limited be- 
cause of the unavailability of direct observations. For the moment, 
the best constraints on esF.Popiii come from the NIRB. 

Since epopiii does not affect the low redshift observations, the 
most stringent limits on it come from the constraints on t^i. In 
fact, matching the WMAP constraint requires that this efficiency 
is confined in the range 0.002 < epopm < 0.03. We would like 
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to stress that, for the assumed value of this constraint on 

epopiii is independent of the constraints obtained from the NIRB 
data. The lower limit on epopiii implies that in order to match 
the WMAP constraints, one requires a star-forming efficiency of 
esF.Popiii = 0.2% if /esc = 1 (probably reasonable for low- 
mass haloes), while in case the escape fraction is much lower, say, 
/esc = 10% (probably true for high-mass haloes), then one needs a 
star-forming efficiency as high as 2%. Also note that the above lim- 
its are obtained for a transition redshift of = 10. The value 
of ztrans is qulte Well constrained to be > 9 from NIRB studies. 
However, for the sake of completeness, we would like to discuss 
the dependence of these limits on ztrans. For higher (lower) values 
of epopiii, one has to take higher (lower) values of Ztrans to obtain 
the same value of Td. Physically, this implies that one can either let 
the PopIII stars form efficiently but survive for a shorter time, or 
let them form inefficiently but survive longer, so that the number of 



free electrons produced is similar. As noted from NIRB studies, the 
value of ztrans wc are using is somewhat a lower limit. This implies 
that the lower limit of 0.002 on epopiii (the case illustrated in Fig. 
4) is quite solid. On the other hand, one can allow for values larger 
that 0.02 if the value of ztrans is allowed to be larger Thus, given 
the present observational constraints, it is slightly difficult to put 
a tight upper limit on epopiii with uncertainties in ^trans. On the 
other hand, if we fix the value of 2;trans from NIRB (or other stud- 
ies in future), our model can constrain the value of epopiii quite 
stringently. 

7.3 Reionization constraints 

In this section we study the constraints on the reionization of both 
hydrogen and Hell with respect to the uncertainties in various free 
parameters of the fiducial model. We have seen that the fiducial 
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model predicts that Hell was first completely reionized at 2: ~ 12. 
However, it turns out that the values of ztrans and epopiii can affect 
the Hell reionization at such high redshifts. Since Hell reionization 
is not complete before z ~ 12, it is obvious that a value of ztrans > 
12 will prohibit a complete reionization of Hell. However, as is 
seen from NIRB studies, it is most likely that the value of ^trans 
may not be much larger than 10. On the other hand, for smaller 
values of epopiii, it is possible that there are not enough photons 
from PopIII stars and the Hell reionization at high redshifts is not 
complete. However, in the extreme case where epopiii = 0.002 
(the lowest allowed value from WMAP constraints), we find that 
QHeiii achieves a maximum value of unity around z ~ 10 (see 
Fig. 4). This implies that Hell must be completely reionized at high 
redshifts as long as the value of epopiii does not violate the WMAP 
constraint. 

Although our fiducial model predicts Hell recombination af- 



ter its first reionization followed by a second reionization driven by 
QSOs ionizing power, a similar double reionization of H does not 
appear to have occurred according to our analysis. The physical in- 
terpretation of Hell recombination at 3.5 < 2: < 7 is found in the 
rather abrupt halt in the PopIII star formation activity at 2 ~ 9; 
as a consequence, the total number of Hell-ionizing photons drops 
considerably and the second Hell reionization has to wait for the 
newly available ones produced by QSOs. This scenario seems to 
strongly imply a double Hell reionization (keeping in mind that the 
first reionization might not be complete for lower values epopiii 
or for higher values of 2trans). Hydrogen, however, behaves differ- 
ently due to the larger availability of photons with energies above 
1 Ryd, which are produced also by galaxies in addition to QSOs. 
Hence the photoionization rate of H remains high enough to keep H 
atoms in the ionized state. A double H reionization could only oc- 
cur if the time gap between the turn-on of PopIII stars and the raise 
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of PopII ones in increased, i.e., if a larger ztrans value is assumed. 
This case has been explored and is shown in Fig. 5, where we have 
fixed ^trans = 11.4. Such model predicts a double H reionization 
around a = 8, but it is at odd with the constraints from the NIRB, 
which require PopIII stars down to redshift of about 9 in order to fit 
the background intensity (and fluctuations) observed in the J band 
(Salvaterra & Ferrara 2003; Magliocchetti, Salvaterra, & Ferrara 
2003). This scenario is a vanilla one for future IGM detection via 
HI 21-cm line, as the redshifted emission could be observed in the 
best range of frequencies and sensitivities of radio telescopes Uke 
LOFAR^ (Ciardi & Madau 2003; lUev et al. 2002). 



7.4 Different radiative feedback models 

In this section, we analyze the effects of varying the strength of 
the radiative feedback we have imposed to the reionization pro- 
cess. As we have already noted, the radiative feedback affects the 
Hell reionization as well as star formation rate at high redshifts. 
Hence it is important to check how the results vary once we change 
our assumptions regarding the feedback mechanism. In our fidu- 
cial model, we have taken the minimum circular velocity of star- 
forming haloes in the ionized regions to be evolving depending on 
the temperature of the region [see equation (14)]. However, a dif- 
ferent feedback prescription is to fix the minimum circular velocity 
of haloes allowed to form stars at a given value, which we take to be 
Vc = 35 km s^\ following (Gnedin 2000; Kitayama et al. 2001). 

We have first explored such case. However, a comparison with 
the fiducial model, has shown hardly any significant differences. 
Therefore we have further increased Vc to 50 km s~^ (a value at the 
upper limit of the physically admissible range derived by the above 
mentioned studies). The corresponding result is shown in Fig. 6. 
While most of the results are unaffected, we find that the strong 
radiative feedback suppresses the growth of Helll regions and de- 
lays the (first) reionization of Hell. As one can see from Panel (b), 
there is substantial decrease in Qneiii when Qhii becomes unity, 
thus delaying the Hell reionization until z w 10. Even this ex- 
treme case, however, does not change the qualitative evolution of 
reionization with respect to the standard one adopted in all other 
cases. Whether this strong radiative feedback has any observational 
consequence is not clear as all the other considered quantities are 
essentially unaffected by the choice of minimum circular velocity. 



7.5 Reduced power on smaller scales 

Finally, we briefly mention some of the other physical processes 
which could reduce the power of density fluctuations on smaller 
scales and hence influence our results. The first such process is 
the suppression of cooling in mini-haloes (i.e., haloes having virial 
temperatures less than 10^ K) due to photodissociation of molecular 
hydrogen. It is also possible that the minihaloes do not contribute 
to reionization due to confinement of HII regions from feedback 
effects (Ricotti, Gnedin, & Shull 2002). This would increase ef- 
fectively the value of Mmiri(z) used in equation (8). Since these 
minihaloes can potentially contribute a substantial number of ion- 
izing photons from PopIII stars at high redshifts, it is obvious that 
one has to use a relatively higher value of epopiii in the case when 
the minihaloes are not forming stars so as to match the WMAP con- 
straints. It turns out that one requires epopiii > 0.007 to match the 

^ http://www.lofar.org 



WMAP limits. This implies that if the high mass haloes have an es- 
cape fraction /eac = 10%, then the required value of star-forming 
efficiency would be as high as 7%, while lower values of the escape 
fraction would require higher star-forming efficiencies. Once such 
a higher value of cpopiii is chosen, we found that there is hardly any 
difference in the reionization history of hydrogen when compared 
to the fiducial model. The reionization of Hell at high redshifts, 
however, occurs much earlier because of less severe feedback from 
ionized HII regions; in fact. Hell reionizes almost simultaneously 
with hydrogen. To understand this, note that the effect of feedback 
is to suppress the photon-production in low mass haloes (which 
were capable of producing photons before the medium was heated 
up) in the ionized (and heated) regions - hence it should be obvi- 
ous that the feedback is more severe when there are more low mass 
haloes forming stars. It is precisely because of this reason that the 
feedback plays a relatively sub-dominant role when the number of 
low mass star-forming haloes is reduced. Feedback effects are in- 
stead more severe in our fiducial model, where the low mass haloes 
are allowed to form stars efficiently. Interestingly, a similar effect 
is produced by smaller values of the fluctuation spectrum index n, 
as a result of the reduced small scale power. 



8 SUMMARY AND CONCLUSIONS 

We have developed a simple formalism to study cosmic reioniza- 
tion and the thermal history of the intergalactic medium (IGM). In 
spite of its simplicity, the formalism implements most of the rel- 
evant physics governing these processes in a self-consistent man- 
ner, including the inhomogeneous IGM density distribution, three 
different classes of ionizing photon sources (massive PopIII stars, 
PopII stars and QSOs), and radiative feedback inhibiting the forma- 
tion of stars in galaxies below a certain circular velocity threshold. 
Such approach allows us to predict the star formation/emissivity 
history of the source, follow the evolution of H and He reionization 
and of the intergalactic gas temperature, along with a number of 
additional predictions involving directly observable quantities. 

By comparing these results with the available experimental 
data we have selected a "fiducial" model. This fiducial model self- 
consistently predicts values matching very well all the available ob- 
servational data, i.e., the redshift evolution of Lyman-Umit absorp- 
tion systems, Gunn-Peterson and electron scattering optical depths, 
and the cosmic star formation history without the need to further 
adjust the free parameters, which are essentially the efficiencies 
of ionizing photon production for PopIII and PopII stars denoted 
by epopiii and epopii, respectively. In principle, a third parameter, 
^trans, the epoch of cosmic transition from PopIII to PopII stars en- 
ters the calculations, but in practice this values is bound to values 
very close to a ^ 9 by the analysis of the Near InfraRed Back- 
ground (NIRB) data. 

The emerging scenario from our analysis can be summarized 
in a few points: 

• Hydrogen reionization must have taken place at z w 15, while 
Hell must have been reionized by z ~ 12, taking into account 
the maximum uncertainty in the value of epopiii (we recall that 
the ionizing photon efficiency is defined as the product of the star 
formation efficiency and the escape fraction, e = esp/esc for each 
Population). At about z = 7 Helll suffered an almost complete 
recombination (QhoIii ~ 0-4) as a result of the extinction of PopIII 
stars (the only contributors of photons with energies above 4 Ryd 
at that redshift), an occurrence required by the interpretation of the 
NIRB. A complete reionization of Helll occurs at = 3.5 and it is 
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driven by hard photons produced by QSOs. Double H reionization 
does not take place due to the larger availability of photons above 
1 Ryd from PopII stars in galaxies and from QSOs, even after all 
PopIII stars disappeared. 

• Following the first reionization, the temperature of the IGM 
corresponding to the mean gas density. To, is boosted to 1.5 x 
10* K; from there it decreases because of adiabatic expansion; 
however, the overall trend is relatively flat. Observations are con- 
sistent with the predicted temperature of the HII regions at 2: > 3.5, 
while they are consistent with the temperature of the Helll regions 
atz < 3.5. This could be interpreted as a signature for the (second) 
reionization of Hell; however, the global (mass-averaged) temper- 
ature rises gradually rather than showing a sudden jump. This al- 
leged jump in the data (Schaye et al. 2000) might be a spurious 
feature induced by fluctuations along the limited number of lines 
of sight used by the absorption line experiments, or it is possible 
that the jump does not have a strong enough statistical significance 
(e.g., see the analysis by Ricotti, Gnedin, & ShuU 2000). 

• PopIII stars produce a first maximum in the star forming ac- 
tivity at z w 15 where psr w 0.003M© yr~^ Mpc~*, followed 
by a drop and a subsequent raise due to the increasing contri- 
bution of PopII stars leading to a less pronounced peak, psF ~ 
O.IMq yr^^ Mpc~^ at z = 4. These results also suggest that only 
a 0.3% of the stars produced by 2: = 2 need to be PopIII stars in 
order to achieve the first reionization. We point out that the fiducial 
model not only correctly predicts at the same time the IGM reion- 
ization and thermal history but it reproduces (with the same free 
parameters) the cosmic star formation history data at 2; < 6. 

In addition to the above features, the data yield the following 
constraints on these free parameters: epopii < 0.01, 0.002 < 
epopiii < 0.03. Varying the efficiencies in these two ranges does 
not affect most of the general scenario emerging from the first two 
points above. We have experimented with maximally strong radia- 
tive feedback finding that the only difference with respect to the 
fiducial case is to suppress the growth of Helll regions at z > 10 
and thus delay the reionization of Hell. Feedback effects impact 
the ionization history of Hell less severely when the minihaloes 
(i.e. haloes with virial temperature T^it < 10* K) at high red- 
shifts are not able to cool and form stars, or when the index of 
the density power spectrum is smaller. However, when the contri- 
bution from the minihaloes is ignored, one requires a higher value 
of star-forming efficiency to match the WMAP constraints. 

Before we concluding we critically discuss our results in the 
light of previous numerical works and semi-analjftical arguments 
in the literature. A few full numerical simulations of cosmic reion- 
ization including radiative transfer have been performed after the 
WMAP results (Ciardi, Ferrara, & White 2003; Ricotti & Ostriker 
2004a; Sokasian et al. 2003; Sokasian et al. 2004; Gnedin 2004). 
These works use a very different set of assumptions concerning the 
IMF of the first stars, feedback effects and recipes for the sub-grid 
physics and radiative transfer schemes. Therefore, a detailed com- 
parison is not possible for all these works. Ciardi, Ferrara, & White 
(2003) concentrated on reproducing the WMAP Thomson optical 
depth in a reionization model dominated by "normal" (i.e. not very 
massive) PopIII stars forming in objects with Tvir > 10* K: mini- 
haloes are thus suppressed (according to our results this should not 
make a sensible difference) but no radiative feedback is included 
for objects just above this threshold. This results in a SFR nearly 10 
times larger than we predict for z > 10. Because of the compen- 
sating effects of IMF and star formation history, hydrogen reion- 
ization is achieved at about the same redshift (14.7) as in our case. 
It is not clear if this model is able to explain the Gunn-Peterson 



opacity as the simulation ends at 2: ~ 8. Gnedin (2004) took the 
opposite approach to the problem, focussing on the fit to the Lya 
mean transmitted flux. The value of esF is fixed by normalizing 
the SFR to its observed value at z = 4, and therefore, as we also 
satisfy this constraint, it should be similar to the one adopted here. 
The production of ionizing photons (which includes our choice of 
the IMF and of the escape fraction) is set by construction to the 
value allowing the best fit to the transmitted flux. Hydrogen reion- 
ization occurs at 2; = 6.1 ± 0.3 in the fiducial model of Gnedin 
(2004) which corresponds to Tci = 0.06 - quite outside the canon- 
ical WMAP value Tei = 0.17 ± 0.04. Although the quoted error 
is somewhat uncertain,* such a low value for Tei seems to be very 
unlikely. For this reason Gnedin suggests that the reionization his- 
tory prior to the hydrogen reionization must have been much more 
complex than the smooth, monotonic behavior obtained from his 
simulations. Our results do not require such a complex history. The 
observed Gunn-Peterson optical depth raise towards 2; = 6 is sim- 
ply caused by the drop of the photoionization rate [see Panel (i) 
of Fig. 1] following the disappearance of PopIII stars, thus caus- 
ing a significant increase of the hydrogen neutral fraction. To be 
more quantitative, we find that the mass-averaged neutral fraction 
climbs from w 2 x 10~* at z = 4 to a maximum of « 5 x 10~* 
at z — 6; following that it decreases until 2: = 9 where it reaches 
the value ~ 4 x 10^''. This behaviour closely tracks (but with op- 
posite derivative) the evolution of the photoionization rate in the 
fiducial model [see Panel (i) of Fig. 1]. It is also consistent with 
the trend, found e.g. by Fan et al. (2002), which explains the good 
agreement with the Gunn-Peterson optical depth of Panel (k). The 
actual values we find might indicate at 2; « 6 a slightly less neutral 
medium than found by Fan et al. (2002) using Lya line - the most 
reliable determination is xni > 5 x 10~* when mass-averaged, 
while xm > 2 x 10~* when volume-averaged. To compare with 
the volume-averaged quantity, we can use our predictions corre- 
sponding to the mean density, which can be shown to yield the 
same answer within about 15%. For the volume-averaged neutral 
fraction, we find a value in the range 1 — 2 x 10~*, which is in 
good agreement with the data. Since Fan et al. (2002) state that the 
volume-averaged fraction should be considered as more reliable, 
we can consider this prediction as a success of our model. 

Two more arguments that have been put forward concerning 
reionization history are worth discussing. The first point concerns 
a higher value of the lower limit (10%) for the neutral fraction de- 
rived by Wyithe & Loeb (2004), based on the size of the ioniz- 
ing radiation influence region around two QSOs at z w 6.3 (see 
also Mesinger, Haiman, & Cen 2004). In brief, the argument used 
is that the size of the ionized region derived from observations is 
smaller than what expected if xm ^ 10^'^, as we argue above. A 
considerable number of uncertainties could jeopardize this conclu- 
sion: (a) the lifetime of QSOs; (b) effect of peculiar velocities; (c) 
radiative transfer (shadowing) effects occurring in the dense envi- 
ronment surrounding high-redshift QSOs (d) determination of the 

* The uncertainty quoted for this number depends on the analysis technique 
employed. Fitting the TE cross power spectrum to ACDM models in which 
all parameters except t^i assume their best fit values based on the TT power 
spectrum, Kogut et al. (2003) obtain a 68% confidence range, 0.13 < t^i < 
0.21 (the one adopted in this paper). Fitting all parameters simultaneously 
to the TT and the TE data, Spergel et al. (2003) obtain 0.095 < t<,i < 0.24. 
Including additional data external to WMAP, these authors were able to 
shrink their confidence interval to 0.11 < Tei < 0.23. Finally, by assuming 
that the observed TT power spectrum is scattered to produce the observed 
TE cross-power spectrum Kogut et al. (2003) infer 0.12 < Tei < 0.20. 
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actual HII region size from the spectrum. For these reasons we feel 
that using this argument as a constrain to reionization models is still 
premature. 

The second point concerns the temperature evolution of the IGM. 
Hui & Haiman (2003), following the original proposal by Theuns 
et al. (2002), noted that as long as the universe is reionized before 
« = 10, and remains highly ionized thereafter, the IGM reaches an 
asymptotic thermal state which is too cold compared to observa- 
tions at 2 = 2 — 4. This indeed applies to our fiducial case, which 
predict a redshift of hydrogen reionization ~ 15. There not only 
we find that the temperature is decreasing relatively slowly due to 
the photoheating provided by He complex reionization history, but 
also that global temperature rises gradually (following the smooth 
QHeiii growth) rather than showing a sudden jump. We have spec- 
ulated that the observed sudden jump in To (i e- the temperature 
of mean density gas) might be a spurious effect due to the limited 
available number lines of sight. This hypothesis could be checked 
in the future by using our results to generate different lines of sight 
to be compared directly with observations. Also, it is possible that 
the statistical significance is much weaker, particularly if the error- 
bars on To are larger than shown here (Ricotti, Gnedin, & ShuU 
2000). 

As a final note, it is useful to remember that our study has not 
included one possible additional contribution to the ionizing back- 
ground due thermal emission from gas shock heated during cos- 
mic structure formation, recently suggested by Miniati et al. (2004). 
Such emission is characterized by a hard spectrum extending well 
beyond 4 Ryd, and according to that study, it is comparable to the 
QSO intensity at redshift > 3. Thermal photons alone could be 
enough to produce and sustain He II reionization already at 2 = 6. 
If this prediction is correct, the partial recombination of Helll seen 
between redshift 3.5 and 8 in our fiducial model might be prevented 
by such radiation. The present results make the test of the He state 
at these intermediate redshifts a crucial benchmark to assess the 
importance of such thermal emission. 
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